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I.   INTRODUCTION 

Computational  fluid  dynamics  has  matured  significantly 
within  the  past  decade  because  of  the  development  of 
increased  computational  capabilities  and  powerful 
computational  techniques.  Current  problems  being  addressed 
include  predicting  flow  separation  over  airfoils  and 
post-stall  flight  characteristics.  These  areas  are  of 
interest  because  studies  [Ref.  1]  indicate  increased  lift 
and  thus  sustained  flight  are  attainable  when  an  airfoil  is 
dynamically  stalled,  that  is,  its  angle  of  attack  is  pitched 
to  a  post  stall  angle  of  attack  rather  than  being  initially 
placed  at  that  high  lift. 

Computational  methods  utilizing  the  full  Navier-Stokes 
(N-S)  equations  are  capable  of  addressing  these  issues,  as 
are  methods  that  include  approximations  to  the  Navier-Stokes 
equations.  One  method,  the  Interactive  Boundary  Layer  (IBL) 
technique,  developed  by  Tuncer  Cebeci  at  Douglas  Aircraft 
Company  and  at  the  California  State  University  [Ref.  2], 
divides  the  flow  over  an  airfoil  into  a  viscous  inner 
boundary  layer  and  an  inviscid  outer  layer.  The 
characteristics  of  the  inner  flow  are  obtained  from  a 
numerical  solution  of  Prandtl's  boundary  layer  equation  and 
the  outer  flow's  characteristics  are  determined  from  Hess 
and  Smith's  panel  method,  and  Fourier  analysis  and  conformal 


mapping.  The  inner  and  outer  layers  are  then  redetermined 
by  an  interaction  model  that  iterates  between  the  two 
regions  and  marches  downstream  until  the  flow  conditions 
have  been  satisfied  at  the  boundary  for  both  regions.  The 
Cebeci  IBL  code  uses  Michel ' s  criterion  to  predict 
transition  from  laminar  to  turbulent  flow  or  transition  may 
be  prescribed.  An  algebraic  (Cebeci-Smith)  turbulence  model 
is  used. 

A  full  Navier-Stokes  code  developed  by  N.L.  Sankar  and 
his  associates  at  the  Georgia  Institute  of  Technology  [Ref. 
3]  uses  an  implicit  finite-difference  procedure  to  solve  the 
2-D  Reynolds-averaged  compressible  Navier-Stokes  equations 
in  strong  conservative  form.  The  time-marching  algorithm 
used  is  an  Alternating  Direction  Implicit  (ADI)  procedure 
developed  by  Beam  and  Warming  [Ref.  4]  and  implemented  by 
Steger  [Ref.  5]  .  The  Sankar  N-S  code  uses  a  body-fitted 
C-grid  system  and  an  algebraic  (Baldwin-Lomax)  turbulence 
model . 

The  Cebeci  IBL  and  the  Sankar  N-S  codes  are  designed  for 
different  purposes.  A  low  Reynolds  number  flow  over  an 
airfoil  tends  to  be  laminar  until  separation.  The  flow  then 
transitions  to  turbulent  flow  and  reattaches  as  turbulent 
flow.  The  Cebeci  IBL  code  models  this  separation  bubble  if 
transition  is  specified  within  the  separation  bubble. 
Velocity  profiles  and  skin  coefficients  are  extremely 
important  in  analyzing  these  low  Reynolds  flows.   Cebeci  has 


developed  codes  for  compressible  oscillating  airfoils, 
however,  this  Cebeci  IBL  code  was  developed  for 
incompressible  steady  state  flow  only  and  thus  does  not 
predict  the  effects  of  unsteady  flow  nor  compressibility. 

The  Sankar  N-S  code  was  developed  to  address  dynamic 
stall  and  its  implication  of  increased  lift.  Therefore,  the 
values  of  interest  to  date  have  been  coefficients  of 
pressure,  lift,  moment,  and  the  effect  of  hysteresis  on 
these  values.  However,  the  Sankar  N-S  code  assumes  the  flow 
is  fully  turbulent,  and  therefore  does  not  account  for 
transition  from  laminar  to  turbulent  flow. 

Neither  dynamic  stall  nor  transition  within  a  separation 
bubble  are  easily  quantified  experimentally.  Transition  is 
a  boundary  layer  phenomenon  and  the  velocity  profiles  and 
skin  frictions  within  the  boundary  layer  must  be  measured  to 
assure  correct  interpretation  of  the  flow  under 
investigation;  surface  pressures  are  not  sufficient  to 
accurately  locate  flow  separation  and  reattachment  [Ref .  6] . 
This  is  a  time  consuming,  expensive  process  prone  to  error. 
Experimental  methods  include  laser  anemometry  and  hot  wire 
probes  [Ref.  7].  Disturbance  of  the  boundary  layer  flow  due 
to  probes  is  undesirable  and  hot  wire  probes  are  normally 
unable  to  determine  flow  direction;  therefore  laser 
anemometry,  although  expensive  and  tedious,  is  increasingly 
being  used. 


Dynamic  stall  is  difficult  to  characterize  due  to  the 
transitory  nature  of  the  phenomenon.  Experimental 
techniques  and  apparatus  include  pressure  transducers,  hot 
wire  probes,  and  laser  doppler  velocimetry  [Ref.  8].  Flow 
visualization  is  also  a  very  effective  tool  for  studying 
both  dynamic  stall  and  separation  bubbles  [Refs.  9,10]. 

A  good  review  of  the  current  state-of-the-art 
computational  and  experimental  aspects  of  aerodynamic  flows 
is  given  in  the  proceedings  of  three  symposia  on  this  topic 
edited  by  T.  Cebeci  [Ref.  11] . 


II.   OBJECTIVES 

The  intent  of  this  study  is  two-fold:  to  become 
familiar  with  computational  fluid  dynamic  methods  and  to 
evaluate  two  codes  to  determine  their  range  of 
applicability. 

Computational  fluid  dynamics  consists  of  various 
mathematical  methods  and  implementation  schemes.  A 
significant  portion  of  analysis  inherent  in  computational 
codes  is  empirical;  therefore,  the  assumptions  used  strongly 
influence  the  results.  It  is  important,  when  attempting  to 
choose  a  computational  code  for  a  specific  purpose,  to  be 
familiar  with  the  significance  of  the  analytical  methods, 
assumptions  made,  and  empirical  models.  Each  code  is 
different  in  these  respects  and  must  be  analyzed 
individually  and  in  detail  to  assure  reliable,  accurate 
results,  especially  when  extending  the  flow  regime  or 
airfoil  to  conditions  whose  features  are  unknown. 

The  two  codes  chosen,  the  Cebeci  IBL  code  and  the  Sankar 
N-S  code,  are  a  good  representation  of  two  powerful, 
accurate  methods  that  differ  widely  in  computational 
approach.  The  Cebeci  IBL  code  has  been  extensively  tested 
in  a  variety  of  steady  state  conditions  [Ref.  12]  and  the 
Sankar  N-S  code  has  compared  well  to  experimental  data  for  a 
pitching  airfoil  [Ref.  13].   However,  the  lack  of  steady 


state  boundary  layer  data  for  the  Sankar  N-S  code  indicated 
that  a  more  in-depth  analysis  of  the  applicability  of  the 
code  was  required.  The  Cebeci  IBL  code  was  used  as  the 
reference  for  the  Sankar  code. 

The  analysis  included  the  following: 

1.  Assess  Ci  for  Reynolds  numbers  of  1.5  and  6  million 

2.  Assess  Cp  and  Cf  for  a  range  of  Reynolds  numbers  (1-15 
million  at  0  degrees  angle  of  attack)  and  angles  of 
attack  (0,  2,  4,  and  6  degrees  for  1.5M  Re  number  and 
0,  4,  8,  and  12  degrees  for  6M  Re  number) 

3.  Assess  velocity  profiles  for  Reynolds  numbers  of  1,  6, 
and  15  million 

4.  Assess  the  influence  of  dissipation  factors  and 
grid  size  on  the  results  of  the  Sankar  N-S  code. 


III.  MATHEMATICAL  FORMULATION 

A.   GOVERNING  EQUATIONS 

Flow  over  an  airfoil  can  be  described  by  the  velocity 

S\  s\  s\  /\ 

q  =  ui  +  vj  +  wk,  the  pressure,  the  density,  and  the 
temperature.  These  six  variables  (u,  v,  w,  p,  p,  and  T)  are 
fully  described  by  the  continuity  equation,  the  equation  of 
state,  p  =  pRT;  the  energy  equation,  6Q  -  6w  =  <$E;  and  the 
three  equations  of  motion.   [Ref.  14] 

The  continuity  equation  states  mass  is  conserved; 
i.e.,  the  flux  of  mass  through  a  cube  per  time  is  equal  to 
the  time  rate  of  change  of  mass.  This  is  shown  in  Figure 
3 . 1  for  the  flow  through  the  faces  perpendicular  to  the  x 
axis.   Mathematically  this  is  expressed  as 


[1CM.  Ax]AyAz  -  [3<EL  Ay]AzAx  -  [UgL  Az]AxAy 


8x     J  J     l  dy 

=  |^  (pAxAyAz)  (3.1) 

Since   the   control   volume   is   fixed,  AxAyAz   is 
independent  of  time;  therefore 


|£  +  iiPuL  +  iiPvL    jMtart. .  0  (3.2) 

St  dx  8y  9z 


or 


p  u  Ay  Ax 


Ay 


Ax 


[pu  i3(P^Ax]AyAz 


-e-  x 


Net  flux  through  face  perpendicular  to  x  axis 
=  _  [_9^Pu)  Ax]  AyAz 

d  X 


Source:   [Ref.  14:p.  106] 

Figure  3 . 1  The  Flux  of  Mass  Through  a  Cube 
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|£+  A    •    pq=  0    .  (3.3) 


Alternately, 


i+p<i+l+i'  =  °  (3-4» 


or 


^  +  p(A-q)    =  0  (3.5) 


The  three  equations  of  motion,  one  for  each  axis  of  the 
Cartesian  coordinate  system,  are  described  by  Newton's 
second  law,  AF  =  A  (ma)  .  The  summation  of  the  x  components 
of   surface   forces   on  the  element  shown   in  Figure   3.2    is 


AF     =  Amax  =    (pAxAyAz)  ^ 


EH:  (3.6) 


with 


do 
AF    =  -  a  AyAz  +  ■  (a     +  -^  Ax)  AyAz 

XX  X  dX 


8t 

t       AxAz  +    (t       +  -s2=.  Ay)AxAz 

yx  yx         3y 


TzxAxAy+    (Tzx  +  -15rAz)AXAy 


3x 

—  A^AxAv    .  (3.7) 


STRESS-STRAIN  RELATIONS 


Figure   3.2      X-components   of  Surface   Forces   on  an  Element 
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Dividing  by  the  volume  of  the  element  yields 


do  3x     3x 

X       yx  ,     ZX   _  n  OU  tr,     o\ 

Tx  +  Ty  +  "5i"  "  p  Dt  (3'8) 


Similarly,  for  the  y  and  z  directions 


do  3x    3t 

_^1  +  _JSY   +   **  =  p  P^  (3.9) 

3y     3x     3z     K  Dt 


and 


do  3x  3t  _ 

z     xz     yz  Dw  ,-,  nr.x 

— 5 —  +  — 5 +  — sr* —  =  P  rr                (3.10) 

3z  3x      3y  Dt 


For  Newtonian  fluids  with  a  single  viscosity  coefficient, 
the  normal  and  tangential  shear  stresses  are  as  follows 
[Ref.  15]: 


ax  =  -p  +  2y  |H  -  |y  (v-q)  (3.11a) 


ay  =  "P  +  2^  f?  "  I1  (V^} 


(3.11b) 


a   =  -p  +  2y  2j  -  |y  (V-q)  (3.11c) 


z    *    '       *    8z    3 


T     =  T     .  V(|5+  |£)  (3. lid) 

yx     xy      3x    3y 


T     =  T     =  y(f£  +  |X)  (3. lie) 

yz     zy      3y    3z 


T    =  T    =  y(|H  +  ||)  (3.11f) 

zx     xy    p  3z    3x 
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Substituting  these  stress  definitions  into  the  equations  of 
motion,  and  assuming  a  constant  viscosity  corresponding  to 
the  mean  temperature  of  the  fluid  ultimately  yields  the 
Navier-Stokes  equations: 


9x     9y     9z 

-  §  +  "A  +  4  ♦  4i  ♦  §  lyiv-ql  =  p  £      0.12b) 

J      9x     9y     9z         J 

2      2      2 

3p  .   r9  w  ,  9  w  ,  9  w,  ,  y  9  rn.^i    „  Dw         ,-j  10_n 

-  9#  +  y[rr  +  7T  +  rr1  +  f  97[v  ^  =  p  de  (3*12c) 

9x     9y     9z 


or  in  vector  format 

-Vp  +  yV2q  +  ^V(V-q)  =  P  f§  +  P(q-V)q         (3.13) 

[Ref.  14] 

B.   REYNOLDS  STRESSES 

The  Navier-Stokes  equations  are  valid  for  laminar  and 
turbulent  flow.  However,  the  complexity  of  turbulence  has 
made  it  impossible  to  relate  the  motion  of  the  fluid  to  the 
boundary  conditions  and  obtain  an  exact  solution. 
Therefore,  the  turbulence  must  currently  be  modeled.  0. 
Reynolds  divided  the  turbulent  flow  into  a  mean  motion  and 
fluctuating,  or  eddying,  motion  as  follows: 
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u  =  u  +  u'  (3.14a) 

v  =  v  +  v'  (3.14b) 

w  =  w  +  w'  (3.14c) 

p  =  p  +  p'  (3.14d) 

p  =  p  +  p'  (3.14e) 

T  =  T  +  T'  (3.14f) 

where  the  barred  terms  are  the  time-average  of  the  component 
and  the  slashed  terms  are  the  fluctuations.  By  definition, 
the  time  averages  of  all  quantities  describing  the 
fluctuations  are  equal  to  zero: 

u'  =  0,    v'  =  0,    w'  =0,    p1  =  0,   p"'  =  0,    T'  =  0 

Rules  for  operating  on  mean  time-averages  are  given 
below.  F  and  g  are  dependent  variables,  and  s  is  the 
independent  variable  x,  y,  z,  or  t. 

I   =  f 

f  +  g  =  ~F+  g 


r  *  ~g  =  f  *  g.       [Ref .  16] 
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The  stresses  caused  by  the  fluctuations  can  be 
determined  using  the  momentum  theorem.  Consider  an  area  dA 
with  dA  *  pu  *  dt  being  the  mass  of  incompressible  fluid 
passing  through  the  element  in  time  dt.  Thus,  the  flux  of 
momentum  in  the  x  direction  is 

dJx  =  dA  *  pu2  *dt;  (3.15a) 

Correspondingly , 

dJv  =  dA  *   p  uv  *  dt  (3.15b) 


and 


dJz  =  dA  *   puw  *  dt.  (3.15c) 

Calculating  the  time  averages  for  the  fluxes  of  momentum 
per  unit  time  yields: 

dJx  =  dA  p    u2  (3.16a) 

dJy  =  dA  p    uv  (3.16b) 

d~J2   =  dA   p   uw.  (3.16c) 
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Utilizing  the  definition  of  turbulent  flow  and  the  previous 
rules  yields 


dJ  =  dA  •  p(u2  +  u'2) 


(3.17a) 


dJ  =  dA  •  p(u*v  +  u'v') 


(3.17b) 


dJ  =  dA  •  p(u*w  +  u'w') 


(3.17c) 


Dividing  these  rates  of  change  of  momentum  by  area  dA,  we 
obtain  stresses.  The  equal  and  opposite  stresses  exerted  on 
the  area  by  the  surroundings  are  a  normal  stress,  -(u2  + 
u'2),  and  two  shearing  stresses,  -(uv  +  u'v1 )  and  -(uw  + 


u'w').   Thus,  the  superposition  of  fluctuations  on  the  mean 
motion  gives  rise  to  three  additional  stresses 


,2 


ax  =  -pu'  ,  t'  =  -pu'v',  x^  =  -pu'w'  . 


(3.18) 


The  total  stress  tensor  due  to  the  turbulent  velocity 
components  of  the  flow  are 


x'  t 
xy    xz 

a'  t' 

y   yz 


yz 


=  "P 


,2 
u' 

u'v' 

u'w' 

,2 

V1 

u'v' 

v'w1 

,2 

w' 

u'w' 

v'w' 

(3.19) 
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The  presence  of  fluctuations  presents  itself  as  an 
apparent  increase  in  stresses  (viscosity) .  These 
additional  stresses  over  the  mean  or  laminar  stresses  are 
termed  apparent  stresses  or  Reynolds  stresses. 

B.   TURBULENCE  MODELING 

The  presence  of  Reynolds  stresses  in  turbulent  flow 
introduces  additional  unknowns  in  the  Navier-Stokes 
equations.  Therefore,  the  Navier-Stokes  equations, 
continuity,  the  perfect  gas  law,  and  the  energy  equation  are 
no  longer  sufficient  to  completely  define  a  solution.  This 
is  known  as  the  closure  problem  and  is  usually  resolved  by 
turbulence  modeling.   [Ref.  17] 

A  common  method  used  is  to  relate  the  turbulent  stress 
to  the  mean  flow  properties  through  empirically  based 
algebraic  formulas.  An  eddy  viscosity,  v^.,  is  defined  in 
the  same  form  as  the  laminar  viscosity.  Previous  models 
related  surface  boundary  conditions  to  points  in  the  fluid 
away  from  the  boundaries  through  wall  functions.  This 
avoided  modeling  the  direct  influence  of  the  eddy  viscosity; 
however,  it  is  only  applicable  in  regions  where  the  Reynolds 
number  is  high  enough  for  viscous  effects  to  be  unimportant 
or  where  universal  wall  functions  are  well  established.  In 
turbulent  boundary  layers  at  low  Reynolds  numbers,  in 
unsteady  or  in  separated  flows,  or  in  three-dimensional 
flows,  the  flow  close  to  the  wall  must  be  described.  [Ref. 
18] 
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A  common  algebraic  turbulence  model  divides  the  flow 
into  an  inner  and  outer  layer.  The  inner  layer  is  defined 
by  a  modified  mixing  length  formula  that  utilizes  some 
damping  function.  The  outer  layer  includes  the  wake  and 
another  damping  function.   [Ref.  17] 

Other  empirical  methods  currently  in  use,  such  as  the  en 
method,  predict  transition.  The  en  method  is  a  stability 
method,  based  on  linear  stability  theory.  It  assumes 
transition  begins  when  a  small  disturbance  is  introduced  at 
or  below  a  critical  Reynolds  number.  The  transition  is 
amplified  by  e9 .  This  method  allows  greater  generality  of 
the  flow,  however  the  formulation  still  relies  on  empirical 
terms.  [Ref.  19] 

The  accuracy  of  turbulence  models  are  limited  by  the 
accuracy  of  the  empirical  constants.  Caution  must  be  taken 
when  using  a  model  under  different  conditions,  i.e.,  a 
different  flight  regime  or  a  radically  different  airfoil. 
The  turbulence  models  mentioned  above  can  be  fairly  simple; 
a  more  complex  model  is  still  not  generally  usable  because 
of  the  computation  costs  involved  and  the  uncertainty  of  the 
constants. 

The  influence  of  turbulence  and  the  transition  from 
laminar  to  turbulent  flow  on  the  airfoil  need  to  be 
understood  and  accurately  modeled  for  a  good  description  of 
the  flow  over  the  airfoil  to  be  detailed.  As  experimental 
methods  continue  to  improve  and  as  computational  methods 
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utilize  the  data,  improvement  in  the  detail  of  the  flow 
field  and  in  flow  prediction  will  follow. 
1.   Cebeci-Smith  Turbulence  Models 

The  turbulence  model  used  in  the  Cebeci  Interactive 
Boundary  Layer  (IBL)  is  a  simple  algebraic  eddy  viscosity 
expression.  Simple  algebraic  models  seem  to  adequately 
predict  turbulent  flow  for  wall  boundary  layer  flows  in 
which  the  Reynolds  shear  stress  and  frequency  do  not  change 
rapidly.  However,  if  the  rate  of  change  of  shear  stress  or 
the  frequency  is  large,  turbulence  models  are  not  currently 
satisfactory.   [Ref.  2] 

The  Cebeci  IBL  code  utilizes  the  algebraic 
eddy-viscosity  formulation  of  Cebeci  and  Smith  [Ref.  20]. 
The  turbulent  eddy  viscosity,  vt/  for  wall  boundary  flows  is 
defined  by  two  separate  formulas;  one  for  the  inner  region, 
based  on  the  Van  Driest  approach,  and  the  other  for  the 
outer  region,  based  on  a  velocity  defect  approach: 


vt  = 


{0.4y[l  -  exp(-y/A)J}2  l^ly^    for  0  <  y  <  yc 


a    /  (ue  -  u)dy  y^  Y         f  or  yc  <  y  £  S 


(3.20) 


where : 


A  =  26v|  (v  |E)      and 
y  max 


1 
Y  =  2  ' 

1  +  5.5(y/6) 
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The  continuity  of  the  eddy  viscosity  defines  yc;  the 
expression  for  the  inner  region  is  used  outward  from  the 
wall  until  it  agrees  with  the  outer  region,  which  is  then 
used.  Ytr  is  an  intermittency  factor  which  allows  for  a 
transition  region  when  progressing  from  laminar  to  turbulent 
flow.   It  is  given  by 

u3    -1.34  x 

v.  =  1  -  exp[ £—*-  R  (x  -  x.  )  /    ^2-J     (3.21) 

**  G    v2  xtr        te   x,  ue 

Ytr  ^ 

where  the  transitional  Re  number,  Re    =  (uex/v^r)  .   The 

xtr 

empirical  constant  Gytr  ^-s  dependent  on  the  Reynolds  number 
of  the  flow.  High  Reynolds  numbers  flows  indicate  Gy^r  = 
12  00,  lower  Reynolds  number  flows  seem  to  be  better  modeled 
by  lower  values  of  Gytr#   [Ref.  12] 

The  parameter  a  in  the  outer  region  is  given  by 
a  =  0.0163/F2*5  where  F  is  the  ratio  of  the  normal  stress 
turbulent  energy  to  the  shear  stress  turbulent  energy 
evaluated  at  the  point  of  maximum  shear  stress.  This  can  be 
expressed  as 


_  1  (u'   -  v'  ) du/dx  )  (3.22) 

|  -  u^V  8u/8y    )(-u'V)max 

The   ratio   of   the   time-averaged   quantities   are 


assumed  to  be  a  function  of  RT  =  Tv/(-u,vl )max  which  can  be 
represented  by 


19 


B  =    -"'  -v'  =  J  (3.23) 

-  u'V     -p-r       )  21^, 

Thus,  the  expression  for  a  is 

a  0-0168 [Ref.    2] 

[1  -  3  0u/8x)/Ou/3y)]2*5 

(3.24) 

Note  that  the  value  of  Y  has  the  effect  of  reducing 
the  eddy  viscosity  away  from  the  airfoil  surface.    This 
turbulence  model  does  not  take  into  account  the  wake  region, 
nor  is  it  validated  for  separated  flow. 
2.   Baldwin-Lomax  Turbulence  Model 

The  Sankar  Navier-Stokes  (N-S)  code  also  uses  an 
algebraic  eddy  viscosity  model,  the  Balwin-Lomax  Turbulence 
Model  [Ref.  21  ]  .  It  is  based  on  the  Cebeci-Smith  two  layer 
model  [Ref.  2  2]  used  in  the  Cebeci  IBL  code  and  may  be 
expressed  as 

2 

!{.4y[l  -  exp(-y/A)]}   |(i)|      for  0  <  y  <  y 
(3.25) 
•0168(1.6)  F„ake  FKLeb(y'     for  ^0^^ 


where 


A  =  26. 
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The  inner  region  is  the  same  mixing  length  formula 
of  the  Cebeci-Smith  model,  simplified.  No  intermittancy 
factor  is  included  (flow  is  calculated  as  wholly  turbulent) , 
A  is  a  constant  rather  than  being  dependent  on  viscosity  and 
velocity  gradients,  and  the  velocity  profile  (du/dy)  is 
replaced  with  the  product  of  vorticity  and  density. 
[Ref.  13] 

The  outer  region  is  based  on  the  wake  function, 
Fwake  and  the  Klebanoff  intermittancy  factor,  FKleb(¥) • 
[Ref.  23] 


F     .      =  min(y         F       ,.25  y         U^.VF       )    .  (3.26) 

wake  vjrmax    max'         -^max    dir    max' 


FKleb(y»   =   [1  +  5"5    (-3^W  ]  (3-27) 


The  quantities  ymax  and  Fmax  are  the  maximum  values  obtained 
from  the  function 


F(y)  =  y  |u|  (1  -  exp  (-y+/A+) }  (3.28) 


which  is  a  form  of  the  mixing  length  formula  used  in  the 
inner  region. 

U^if  is  the  difference  between  the  maximum  and 
minimum  velocity  of  the  velocity  profile.  Fwake^s  similar 
to  the  y  of  the  Cebeci-Smith  turbulence  model  and  thus  also 
reduces  eddy  viscosity  away  from  the  airfoil  surface. 
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This  model  has  been  used  in  separated  flows  and  in 
the  wake,  however  its  validity  is  not  assured  in  these  regions. 
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IV.   THE  BOUNDARY  LAYER  EQUATION 

Prandtl  clarified  the  influence  of  viscosity  in  high 
Reynolds  flows  by  simplifying  the  Navier-Stokes  equations  to 
yield  approximate  solutions.  He  divided  the  air  flow  over  a 
body  into  two  regions: 

1)  The  region  near  the  surface  where  viscous  forces 
dominate . 

2)  The  rest  of  the  flow  where  inertia  forces  dominate; 
this  region  may  be  considered  frictionless  and 
potential. 

Consider  a  2-D  incompressible  flow  over  a  body.   Most  of 

the  flow  is  moving  at  free  stream  velocity.   However,  at  the 

surface  the  velocity  is  zero,  increasing  to  free  stream  at 

some  distance  from  the  surface  as  shown  in  Figure  4.1.   In 

this  first  region,  called  the  boundary  layer,  the  velocity 

gradient  normal  to  the  wall,  3u/8y,  is  very  large,  as  is  the 

shearing  stress, 


T  =  U   ^  (4.D 


The  two  regions  are  not  distinct,  but  are  usually  divided  at 
the  streamline  where  the  velocity  reaches  99%  of  the  free 
stream  velocity. 

Simplification  of  the  Navier-Stokes  equations  will  be 
accomplished  by  doing  an  order  of  magnitude  analysis  of  each 
term.   The  following  assumptions  are  made: 

23 


1         ° 

I 

Pressure  decrease 


Pressure  increase 


Figure  4.1  Velocity  Profile 
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1)  A  flat  wall  coinciding  with  the  x-direction,  with  the 
y  axis  perpendicular. 

2)  6  «  L 
5/L  <<  1 

where  6  is  boundary  layer  thickness  and  L  is  a  linear 
dimension  of  the  body  so  selected  to  ensure  3u/ 3x  =  1 
under  the  region  in  consideration. 

3)  The  Reynolds  number,  R  =  ULp/  y  =  —  is  large.   [Ref. 
16]  V 

The    Navier-Stokes    equations    are    rewritten    in 
dimensionless  form: 

-  Velocities  are  with  respect  to  the  free  stream  velocity, 
U. 

-  Linear  dimensions  are  with  respect  to  a  characteristic 
length,  L. 

-  Pressure  is  divided  by  pU2. 

Retaining  the  same  symbols  for  the  dimensionless  quantities 
as  for  their  dimensional  counterparts,  and  writing  the  order 
of  magnitude  under  each  term  yields  for  continuity: 


|H+|^=0  (4.2) 

3y   3y 


and   for  the  Navier-Stokes  equations 


2  2 

,.  3u   ,        3u  3p   ,    1,3  u   .    8  u,  (A    -,-,x 

dirx      u3x  +  v3?=-#+R(^2  +  ^2)  (4>3a) 

1  2  1 

1     1       6     =■  6Z     1        -4 

6  ^ 

2  2 

,.  3v   ,        3v  3p  ,    l/3v       3v,  /4    2h) 

d~y     u^  +  v3y-="3y     R^?^  (        } 

1661  6       625         j 
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At  the  outer  edge  of  the  boundary  layer  u  equals  the 
steady  flow  U(x) .  The  viscous  terms  no  longer  dominate  and 
thus,  for  the  outer  flow 


u|U=_19p    (dimensional)  (4.4) 


or  rewritten  in  the  form  of  Bernoulli's  equation 

p  +  1/2  p  IT  =  constant  -  (4.5) 

Returning  again  to  dimensional  quantities,  the 
simplified  Navier-Stokes  equations,  known  as  Prandtl's 
Boundary-Layer  Equations  may  be  written: 


S  +  fJ-0  (4.6a) 

2 

3u  ,   8u     1  dp  ,  ,^3u  ,,  g-v,\ 

U  t? V   V  r—  = j   +  V" t"  (4.6b) 

dx  8y     p  dx     a  2 

with  the  boundary  conditions  y=0:  u=0  v=0;y=  inf: 
u  =  U(x) .  Also,  a  velocity  profile  at  the  initial  section, 
x  =  x0,  must  be  prescribed.   [Ref.  16] 
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V.   THE  SANKAR  NAVIER-STOKES  METHOD 

The  Navier-Stokes  Code  utilizes  the  governing  equations 
in  conservative  form,  a  body  fitted  coordinate  system,  and 
an  Alternating  Direction  Implicit  (ADI)  procedure  [Ref.  23]. 
The  governing  equations  in  conservative  form  have  the 
coefficients  of  the  derivative  terms  constant.  The 
conservative  form  ensures  the  conservation  of  the  physical 
properties. 

A.   GRID  GENERATION 

The  requirements  of  a  grid  in  the  physical  plane  and  in 
the  computational  plane  are  conflicting — therefore  a  grid 
transformation  is  advantageous.  For  ease  in  computation, 
equal  spacing  of  grid  points  is  desirable;  however,  the 
physical  grid  needs  to  be  clustered  so  that  the  boundary 
layer  and  sharply  curving  surfaces  such  as  the  leading  edge 
contain  enough  points  so  as  to  be  adequately  defined.  The 
boundary  conditions  must  be  accurate  and  should  be  contained 
on  rectangular  surfaces  in  the  computational  plane.  Also, 
the  grid  should  be  smooth  with  few  discontinuities. 

The  present  code  uses  an  algebraic  C-grid  which 
generates  a  sheared  parabolic  coordinate  system  [Ref.  23], 
first  proposed  by  Jameson  [Ref.  24]. 
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First,  two  points,  T  and  N,  are  defined  on  the  desired 
airfoil  in  the  X-Z  plane  as  shown  in  Figure  5.1.  These 
points  are  defined  by  the  complex  values 

zT  =  Xiji  +   iyij  (5.1a) 

zN  =  xN  +   iyN  (5.1b) 

respectively.  T  is  located  at  the  trailing  edge  of  the 
airfoil  and  N  is  located  half  way  between  the  leading  edge 
and  the  center  of  curvature  of  the  leading  edge. 

Next,  the  airfoil  in  the  physical  plane  is  transformed 
as  shown  in  Figure  5.2.  A  trailing  edge  vortex  sheet  shape 
is  assumed  to  leave  the  trailing  edge  smoothly  by  running 
tangent  to  the  mean  camber  line  at  the  trailing  edge.  The 
airfoil  and  wake  are  then  mapped  onto  the  plane  by  using  the 
following  transformation, 

X,    =  /  z  -  zN  .  (5.2) 

The  NACA  0012  airfoil  shape  transforms  to  that  shown  in 
Figure  5.3.  Cubic  interpolation  defines  additional  points 
to  smooth  the  surface. 

The  far  field  boundaries  are  mapped  in  the  physical  and 
computational  planes  as  shown  in  Figure  5.4. 
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Figure  5.1   Defined  Points  on  the  Airfoil  Surface 
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Source:    (Ref.  13:p.  19) 

Figure  5.2   Symmetric  Airfoil  in  the  Physical  Plane 
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Source:    (Ref.  13:p.  19) 


Figure  5 . 3   Symmetric  Airfoil  Unwrapped  in 
the  Intermediate  Plane 
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Figure   5.4      Far  Field  Boundaries 
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A  sheared  Cartesian  coordinate  system  is  then 
constructed  in  the  x,  plane.  It  consists  of  straight 
vertical  lines  which  then  contain  specified  clustered 
spacings  defined  by  a  stretching  function.  This  allows  the 
grid  size  to  increase  normal  to  the  surface. 

Lastly,  the  grid  is  mapped  back  to  the  physical  plane. 
The  points  on  the  computational  plane,   £  ,  are  given  by 

C  =  U  in  (5.3) 

and  on  the  physical  plane  by 

x  -  xn  =  &   +    n2  (5.4a) 

y  -  yn  =  2£n  (5.4b) 

[Refs.  13,23]. 

The  present  method  uses  a  grid  containing  161  points  in 
the  E,  direction  and  41  points  in  the  n  direction.  The  final 
grid  in  the  physical  plane  is  shown  in  Figure  5.5  and  a 
detail  of  the  grid  around  a  NACA  0012  airfoil  is  presented 
in  Figure  5.6. 

B.   INITIAL  CONDITIONS  AND  BOUNDARY  CONDITIONS 

The  initial  conditions  for  viscous  flows  are  the  free 
stream  conditions.  Viscous  dissipation  inherent  in  the 
equations  will  minimize  errors  in  this  approximation  after  a 
sufficient  time.  Inviscid  flows  require  the  proper 
combination  of  "artificial"  dissipation  and  boundary 
conditions  for  the  correct  result. 
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Figure  5.5   Final  Grid  in  the  Physical  Plane 
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Figure  5.6   Detail  of  Grid  Around  the  Airfoil 
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Boundary  conditions  must  be  defined  on  the  airfoil,  at 
the  far  field,  and  in  the  wake.  On  the  solid  surface,  two 
conditions  must  be  satisfied;  no  penetration  of  the  surface 
and  the  solid  and  fluid  have  the  same  velocity  at  their 
boundary  (the  no  slip  condition) .  Adiabatic  flow  is 
assumed. 

In  the  far  field,  the  flow  is  represented  by  the  linear 
small  disturbance  equation 

(1  -  M^^XX  +  $yy  =  0  (5.5) 

where  $  is  the  perturbation  potential  and  x  and  y  are  the 
physical  plane  coordinates.  This  model  is  used  instead  of 
specifying  flow  conditions  at  infinity  to  compensate  for  the 
loss  of  lift  experienced  when  the  boundary  is  not  placed  far 
enough  away  from  the  solid  surface  [Ref.  23]. 

In  the  wake,  the  grid  procedure  produces  a  "cut"  along 
the  coordinate  line  from  the  trailing  edge  to  the  far  field 
boundary.  Here,  the  flow  properties  are  averaged  from  above 
and  below.   [Refs.  3,13] 

C.   NUMERICAL  FORMULATION 

The  coupled  and  non-linear  governing  equations  are 
solved  using  an  alternating  direction  implicit  time  marching 
procedure  similar  to  that  developed  by  Beam  and  Warming 
[Ref.  4]  and  as  used  by  Steger  [Ref.  5].  The  governing 
equations  are  assumed  to  have  a  solution  at  some  time  tn  and 
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the  solution  is  then  advanced  to  some  time  tn+i  using  Euler 
implicit  time  differencing.  The  equation  is  then  linearized 
using  Taylor  expansion  and  reduced  to  a  series  of 
one-dimensional  problems  using  the  Beam  and  Warming 
approximate  factorization  method  [Ref .  5] . 

The  procedure  is  not  fully  implicit  since  the  viscous 
terms  are  lagged  by  one  time  step.  Artificial  dissipation 
is  included  in  both  second  and  fourth  order  terms  to  obtain 
accurate  surface  pressures.  The  ADI  procedure  is  limited  in 
time  step  because  of  accuracy  but  the  explicit  boundary 
conditions  further  limit  the  step  size  due  to  stability 
considerations . 

1.   Governing  Equations 

The  governing  equation  in  the  computational  plane  is 
given  as 

3Tq  +   3^  E  +  8  F  =  Re_1(8^R  +  dr)   S)  (5.6) 


where 


q  =  q/J  (5.7) 


E  =    Ut<3  +    ?xE   +    CyF)/J  (5.8) 

F   =    (ntq   +    nxE"    +    nyF)/J  (5.9) 
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R  =    (£XR  +   SyS/J  (5.10) 


-*■  -> 


s  =   (nxR  +   nyS)/J  (5.11) 

txx   =    (X+2y)  Uxu£+%un)    +   XCCyV^+riyV^  (5ol2a) 

Txy  =    yCUyU^  +  nyi^)    +    (  ^v^+t^v^)]  (5.12b) 

Tyy   =    (X+2y)  (CyV^+TyV^)     +    X  (£  XU^+  T^U^)  (5.12c) 

R4   =  uxxx  +  VT^  +  kPr"1(t-l)  (^a2+nx3na2)  (5.13a) 

S4   =  uxxy  +  vxyy  +  kPr" 1  ( Y - 1 )  ( £y3 £  a 2 -HI  y  d^ a 2 )  (5.13b) 

J  is  the  transformation  Jacobian: 

J  =  Sxny-£ynx  =   (x^y^x^y^)-1  (5.14) 

Also, 


Sx  =  JYn  (5.15a) 


£y  =   -Jxn  (5.15b) 

nX  =   "Jyr  (5.15c) 
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n  y  =  Jx^  (5.15d) 

St  =  "xT^x-yT^y  (5.15e) 

nt  =  'V^x-Y^y  (5.15f) 

Details  of  the  derivation  of  the  governing  equations  in  the 
physical  plane  and  transformation  of  the  governing  equations 
to  the  computational  plane  are  contained  in  References  13 
and  2  3 . 

2 .   Time  Differencing  and  Linearization 

The  governing  equations  cannot  be  solved  directly  in 
the  form  of  Equation  (5.6)  because  of  their  non-linearity. 
This  is  overcome  by  writing  the  flow  quantities  p,  pu,  pv, 
and  e  at  their  new  time  level  as  their  value  at  the  known 
time  level  and  their  increment,  i.e., 


n+1     n  ,  A.  n+1  /c  t  *-\ 

p     =  p  +  Ap  (5. 16) 


The  variable  q  can  therefore  be  written  at  the  new 
time  level  tn+1  as 

qn+l  =   qn  +  /\qn+l  (5.17a) 


or 
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Aqn+1  =  qn+l  _  qn  #  (5.17b) 

A  Taylor  series  expansion  of  qn  backward  about  qn+1  yields 

qn  =  qn+l  _  At6tqn+1  +  0(At2)  (5.18a) 

qn+l  _  qn  m   At6tqn+1  +  0(At2)  =  Aqn+1  (5.18b) 

Using  the  same  procedure  with  central  differencing  for  the 
spatial  derivatives  and  substituting  into  the  above  equation 
yields 


Aqn+1  =  -At(6rEn+1  +  6  Fn+1)  +  ARe~1(5rRn+1 

+  £  Sn+1)  +  0(At2)_         (5.19) 


where  ^r  and  6^  are  second  order  accurate  difference 
operators . 

Backward  differencing  is  first  order  accurate  and 
central  differencing  is  second  order  accurate,  therefore  the 
Equation  (5.19)  is  first  order  accurate  in  time  and  second 
order  accurate  in  space. 

The  governing  equation,  Equation  (5.6),  is  still 
non-linear.  So  that  it  can  be  solved  directly  rather  than 
iteratively,  it  is  linearized.  First  E  and  F  are  rewritten 
using  a  local  Taylor  expansion  at  qn: 
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En+1  =  En  +  [3qE-|nAqn+l  +  Q(At2)  (5.20a) 

Fn+1  m   Fn  +  [9qF]nAqn+1  +  o(At2)  (5.20b) 

where  [8qE]  and  [ aqF]  are  Jacobian  matrices.  [Ref.  23] 

Substituting  Equation   (5.20)   into  the  time  differenced 
governing  equation,  Equation  (5.19),  yields 


([I]  +  At6  [  aqE]n  +  At6  [8qF]n)<Aq}n+1 

£  n 


=   -At(5rEn+6    Fn)    +    AtRe_1(6rRn+1+6   Sn+1)  (5.21) 


This  can  be  expressed  as 


([I]+At6  [A]  +  At6n[B]){Aq}n+1  =  {R}n  (5.22) 


where : 


[A]  =  [3qE]n  (5.23) 

[B]  =  OqF]n  (5.24) 

{R}n  =  -At(6-En+6  Fn)  +  AtRe_1(6  Rn+l+e^s11"1"1)     (5.25) 

The  governing  equations   for  the  unknown  vector 
{Aq}n+1  are  now  linear  and  may  be  solved  numerically. 
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3.   The  Alternating  Direction  Implicit  Procedure 

An  approximation  factorization  method,  developed  by 
Beam  and  Warming  [Ref.  4]  is  now  used  to  solve  the  governing 
equation  Equation  (5.22),  which  has  been  linearized  in 
[Aq]n+1.  Although  Equation  (5.22)  could  be  solved  directly, 
Beam  and  Warming's  method  reduces  the  large  (and  thus 
costly)  matrix  into  a  sequence  of  one-dimensional  problems. 
The  governing  equation  is  approximated  as 


([I]    +   At5    [A])([I]    +  AtyBmAq}1^1   =    {R}n  (5.26) 


which  is  then  rewritten  as  two  equations 


([I]  +  At5c[A]){Aq*}  =  <R}n  (5.27a) 


([I]  +  At 6  [B]){Aq}n+1  =  {Aq*}  (5.27b) 


Since  6-  and  6  are  central  difference  operators,  Equations 
(5.27)  are  systems  of  block  tridiagonal  matrix  equations 
composed  of  4x4  submatrices.  These  can  be  solved  by  one  of 
several  methods,  such  as  LU  decomposition.  Once  (Aq}n+1  is 
obtained,  {q}n+1  can  also  be  obtained  from 

{q}n+l  =  {q}n  +  {Aq}n+l  (5.28) 
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Boundary  conditions  at  {Aq}  as  well  as  the  defined  vector 
{Aq*}  must  be  defined. 

4.   Stability  and  Accuracy  Considerations 

The  ADI  approach  is  implicit  with  explicit  boundary 
conditions  and  viscous  terms  that  are  lagged  by  one  time 
step.  Implicit  numerical  methods  theoretically  are 
unconditionally  stable  with  the  size  of  the  time  step 
limited  by  accuracy  rather  than  stability.  The  time  step 
stability  limit  imposed  by  the  explicit  boundary  conditions 
must  therefore  be  less  than  the  accuracy  time  step  limit. 
[Ref.  4] 

A  linear  stability  analysis  for  an  explicit 
procedure  is  performed  [Ref.  23]  to  determine  the  maximum 
time  step  possible,  then  a  very  conservative  estimate  is 
used.  When  calculating  steady-state  flow,  the  convergence 
of  the  calculations  can  be  improved  by  introducing  a 
variable  time  step.  This  allows  small  time  steps  where  the 
grid  must  be  highly  clustered,  such  as  the  boundary  layer, 
in  regions  of  shocks,  and  near  stagnation  points,  and  large 
time  steps  elsewhere  where  the  grid  is  larger.   [Ref.  23] 

Viscosity  slows  a  flow  down  by  dissipating  energy. 
Mathematically,  this  results  in  a  reduction  of  flow  field 
gradients.  The  mathematical  formulations  used  to  include 
viscosity  in  a  flow  can  therefore  also  be  applied  to 
suppress  errors  inherently  generated  in  certain  methods. 
[Ref.  25] 
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The  capability  of  including  this  "artificial 
viscosity"  as  used  by  Murman  and  Cole  [Ref.  26]  is 
implemented  through  backward  differencing  in  a  Taylor  series 
expansion.  The  truncation  term  mimics  viscosity  and  is 
known  as  artificial  viscosity.  When  the  method  includes 
this  dissipative  term  it  is  known  as  implicit  artificial 
viscosity.  Often  however,  more  dissipation  is  required  for 
convergence  or  stability,  or  because  it  is  advantageous  to 
apply  selective  dissipation.  Then,  explicit  artificial 
viscosity  or  explicit  dissipation  is  included  in  the 
numerical  formulation. 

Central  differencing  in  a  Taylor  series  expansion 
decouples  the  even  from  the  odd  terms,  causing  high 
frequency  errors.  The  spatial  derivatives  are  formulated 
through  central  differencing  and  thus  contain  these  errors 
which  influence  the  solution  accuracy  when  using  large  time 
steps.  The  time  derivatives  are  formulated  through  backward 
differencing;  their  artificial  viscosity  terms  and  the 
viscous  terms  suppress  this  error  somewhat.  To  further 
dissipate  the  high  frequency  errors,  especially  in  high 
Reynolds  number  and  inviscid  flows,  both  fourth  order  and 
second  order  artificial  dissipation  is  explicitly  applied  to 
the  right  hand  side  of  the  descretized  governing  equations 
as  was  done  by  Jameson  [Ref.  24].   [Ref.  23] 

The  fourth  order  terms  alone  lead  to  overshoots  in 
the  vicinity  of  shocks.   The  second  order  dissipative  terms 
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correct  this  but  tend  to  smear  results  in  regions  such  as 
the  leading  edge.  This  is  resolved  by  implementing  second 
order  dissipation  only  in  regions  of  high  pressure 
gradients.  Including  artificial  dissipation  in  the  right 
hand  side  of  Equation  (5.26)  also  corrects  for  incorrect 
initial  conditions  after  a  sufficient  time  step. 

To  allow  the  viscous  terms  to  be  explicitly 
modulated  and  to  remove  any  explicit  stability  limits, 
artificial  dissipation  is  also  implicitly  included  in  the 
left  hand  side  of  Equation  (5.26).  The  inclusion  of 
explicit  and  implicit  artificial  dissipation  yields  [Ref. 
23] 


([I]  +  At5^  [A)-^"1*  ^T)  ([I]+At6n[B]-?iJ-16nTiJ){Aq} 


n+1 


=  <R}n+1  -  (Di  -  D-j-i)  -  J-16nnnn(Jq)         (5.29) 


where  is  a  function  of  the  maximum  pressure  gradient  and  D 
and  D^_!  are  either  second  or  fourth  order  dissipation.  The 
artificial  dissipation  terms  do  not  affect  the  accuracy  of 
the  formulation  since  all  terms  are  of  the  same  order  as  or 
smaller  than  the  truncation  errors  associated  with  the 
spatial  and  time  difference  formulas. 
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5.   Application  of  Boundary  Conditions 

All  of  the  boundary  conditions  are  explicitly 
applied  after  the  ADI  sweeps  at  each  time  step.  The 
explicit  method  was  used  because  of  its  ease  of 
implementation  in  the  code  even  though  implicit  boundary 
conditions  are  more  desirable  because  of  accuracy  and 
stability  considerations.   [Ref.  23] 

On  the  airfoil  surface  the  no-slip  condition  and  the 
assumption  of  no  flow  through  the  surface  correspond  to 

8p/8n  =  0  3p/3n  =  0  (5.30) 

Adiabatic  flow  is  assumed.  A  two  point  extrapolation  of  the 
above  yield  p  and  p  to  be  [Ref.  14] 

pi,l  =  (4Pi,2  "  pi,3)/3  (5.31a) 

pi,l  =  (4Pi,2  -  pi,3)/3  (5.31b) 

Internal  energy  and  the  coefficient  of  pressure  may  be 
obtained  from  the  calculation  of  p  and  p.  The  incremental 
quantities  {A<3*}  and  {A  q}n+1  are  also  assumed  zero  on  the 
solid  surface  and  are  solved  in  a  similar  manner  as  for  p 
and  p  .  The  far  field  boundary  conditions  are  assumed  to  be 
free  stream  plus  the  disturbances  caused  by  the  far  field 


46 


not  being  placed  at  infinity  as  described  in  Section  V.B. 
Internal  energy  is  evaluated  from 


e  =  -§_  +  .5p(u2  +  v2)  (5.32) 


The  downstream  boundary  conditions  are  extrapolated  from  the 
adjacent  interior  points.  At  the  boundary,  pressure  is 
taken  to  be  freestream. 

The  cut  inherent  in  the  grid  system  divides  the  flow 
above  and  below  the  line  emanating  from  the  trailing  edge. 
The  flow  quantities  on  the  cut  are  obtained  by  averaging  the 
values  of  the  interior  points  above  and  below  the  cut.  This 
is  acceptable  because  of  the  denseness  of  the  grid  in  this 
region.   [Ref.  23] 

D.   USE  OF  SANKAR'S  N-S  CODE 

Reference  13  details  the  use  of  Sankar's  N-S  code  for 
both  steady  and  unsteady  flows.  The  code  is  submitted  to 
the  NASA  X-MP  Cray  via  Job  Control  Cards.  JCL  options  are 
selected  so  to  access  files  where  output  data  is  stored  or 
to  send  data  to  a  specific  directory.  Job  Cards  also 
provide  account  and  time  limit  information. 

The  main  program  contains  all  of  the  "Write"  statements, 
which  may  or  may  not  be  required.  Output  information 
includes  input  data,  airfoil  coordinates  in  the  physical  and 
computational  planes,  grid  information,  residuals,  pressure 
and  skin  friction  profiles,  coefficients  of  lift,  drag  and 
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moment ,  and  velocity  profile  information.  The  time  step 
intervals  at  which  these  are  printed  is  specified  within  the 
program . 

The  majority  of  the  program  inputs  are  located  in  data 
cards.  These  inputs  and  their  definitions  are  decribed  in 
Table  5.1.  See  Reference  13  for  a  detailed  program 
description.  The  values  used  when  implementing  the  Sankar 
N-S  code  are  given  in  Table  5.2. 
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TABLE  5.1 

INPUT  PARAMETERS  FOR  THE  SANKAR  N-S  CODE 

INPUT  DESCRIPTION 

IMAX  Number  of  x  coordinate  locations 

KMAX  Number  of  y  coordinate  locations 

DT  Size  of  time  step 

WW  Explicit  artificial  viscosity  term 

ALFA  Mean  angle  of  attack 

ALFA1  Amplitude  of  oscillation 

ALFAI  Angle  below  which  modified  turbulence 

model  is  used 


REDFRE 

AMINF 

FNSTP 

REYREF 

DNMIN 

TSTART 

FULOUT 

RSTRT 

PITCH 

PLUNGE 

FNU 

FNL 

ZSYN 

XP 

YP 


Reduced  frequency 

Mach  number 

Number  of  time  steps 

Reynolds  number  in  millions 

Distance  of  first  point  off  of  the  wall 

Time  calculation  flag 

Plotting  file  flag 

If  set,  stored  values  are  read  to 
continue  iteration 

Set  for  dynamic  stall,  indicates  change 
in  AOA 

Set  for  up  and  down  motion  of  airfoil 

Number  of  upper  airfoil  coordinates 

Number  of  lower  airfoil  coordinates 

Flag  for  symmetric  airfoil 

X  airfoil  coordinates 

Y  airfoil  coordinates 
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TABLE  5.2 

VALUES 

USED  FOR  SANKAR  STEADY  N-S  CODE 

INPUT 

VALUE 

IMAX 

157 

KMAX 

41 

DT 

1.00 

WW 

10-12 

ALFA 

0-12 

ALFA1 

0 

ALFAI 

0 

REDFRE 

0.0 

AMINF 

.12,  .30 

FNSTP 

2000-8000 

REYREF 

1-15 

DNMIN 

.000005-. 0001 

TSTART 

0 

RSTRT 

0 

PITCH 

FALSE 

PLUNGE 

FALSE 

FNU 

33 

FNL 

33 

ZSYN 

0 

XP 

X  airfoil  coordinat 

0012  airfoil 


YP 


Y  airfoil  coordinates  for  NACA 
0012  airfoil 
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VI.   THE  CEBECI  INTERACTIVE  BOUNDARY  LAYER  METHOD 

The  Cebeci  Interactive  Boundary  Layer  (IBL)  code  uses  a 
two  layer  approach;  an  outer  inviscid  layer  and  a  viscous 
inner  boundary  layer.  In  this  it  follows  Prandtl's  boundary 
layer  theory  by  assuming  the  inner  viscous  flow  is  the  only 
region  where  viscous  effects  are  important.  The  remaining 
outer  region  is  dominated  by  inertia  terms  and  can  be 
assumed  inviscid.  The  inviscid  outer  flow  characteristics 
are  determined  through  Hess  and  Smith's  panel  method.  The 
flow  is  assumed  to  have  a  vortex  and  source  distribution 
such  that  it  gives  correct  circulation  and  velocity  over  the 
airfoil.  The  airfoil  surface  is  replaced  by  a  distribution 
of  panels  that  satisfies  the  Kutta  condition. 

The  inner  viscous  flow  utilizes  the  boundary  layer 
method  to  determine  the  flow  characteristics.  A  direct 
boundary  layer  method  is  used  in  regions  of  large  viscous 
stresses  such  as  near  the  leading  edge.  This  method 
prescribes  an  external  velocity  and  requires  the  no  slip 
condition  to  be  satisfied.  In  the  rest  of  the  flow,  an 
interactive  boundary  layer  method  is  used.  Here,  the  edge 
boundary  conditions  prescribe  a  combination  of  displacement 
thickness  and  external  velocity.  An  iterative  technique  is 
used  to  solve  for  this  flow.   [Ref.  27] 
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The  inner  and  outer  flow  are  coupled  through  an 
interaction  model.  The  boundary  location  and  velocity  are 
unknown  and  are  solved  simultaneously  through  an  iteration 
between  the  inner  and  outer  flow  equations. 

Cebeci's  IBL  code  calculates  where  transition  from 
laminar  to  turbulent  flow  occurs  and  incorporates  a 
smoothing  function.  It  also  predicts  separation  and  allows 
laminar  separation  and  then  transition  to  turbulent  flow  and 
subsequent  reattachment  of  the  turbulent  flow. 

A.   VISCOUS  INNER  FLOW 

1.   Direct  Boundary  Layer  Method 

The  direct  boundary  layer  method  can  be  used  in 
regions  of  the  boundary  layer  where  the  viscous  effects  have 
not  yet  strongly  influenced  the  flow.  This  usually  implies 
the  stagnation  point  and  the  airfoil  leading  edge.  The 
advantage  of  minimal  viscosity  influence  is  the  capability 
of  defining  a  stream  function,  ty,  that  satisfies  the 
continuity  equation 

u  =  di>/dy  and    v  =  -  3^/  3x  (6.1) 

The  momentum  equation  is  subjected  to  the 
Falknar-Skan  transformation  [Ref.  2] 


n  -  Vvl  y  (6*2a) 
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f  (x,n)  =        iMx,y)  (6.2b) 

The  normal  coordinate  y  and  the  stream  function  \\>  are  scaled 
with  respect  to  the  external  velocity  for  convenience  and 
accuracy.  The  boundary  conditions  of  no  flow  penetration 
through  the  wall  and  no  slip  condition  at  the  wall  are  also 
transformed.  The  resulting  momentum  equation  and  boundary 
conditions  are  given  as 

jbf..).  +  l(m+l)ff"  +  m[l-(f')2]  =  x(f  ||!--  f"  ||)         (6.3) 

where : 

b  =  1  +  v^v 

m  =  (x/ue)  (du^/dx) 

and 

n  =  0:    f'(x,0)  =  0,     f(x,0)  =  0 

n  =  ne  :   f »  (xfne)  =  l 

Primes  denote  differentiation  with  respect  to  ri.  This  is  a 
third  order  partial  differential  equation  and  cannot  be 
solved  directly.  Therefore,  the  box  method  developed  by 
Keller  [Ref.  28]  is  used.   [Refs.  12,30] 
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The  box  method  reduces  the  governing  equations  to  a 
first  order  system  through  introduction  of  two  dependent 
variables.  The  flow  properties  are  then  evaluated  only 
discretely  by  defining  the  solution  domain  as  a  rectangular 
mesh.  Instead  of  solving  continuous  functions,  all 
parameters  are  approximated  in  terms  of  nodal  values  and 
their  location  on  the  mesh.  The  domain  of  dependence  is 
substantially  reduced  and  the  overall  solution  scheme 
simplified  by  solving  for  the  nodal  values  through  central 
differencing . 

The  resulting  nonlinear  system  is  solved  by  Newton's 
method.  This  iterative  procedure  linearizes  the  variables 
at  location  i  by  rewriting  the  value  at  i  as  a  sum  of  the 
value  at  location  i-1  plus  some  incremental  value. 
Substituting  into  the  governing  system  of  equations  results 
in  a  linear  system  of  the  unknown  incremental  values  which 
are  repeatedly  solved  until  they  are  small  enough  to  be 
neglected. 

2.   Interactive  Boundary  Layer  Method 

Most  of  the  airfoil  is  influenced  by  viscosity,  thus 
the  direct  boundary  layer  method  can  no  longer  be  used. 
Instead,  an  interactive  boundary  layer  method  is  used.  It 
is  effective  even  in  regions  of  rapid  flow  acceleration, 
boundary  layer  separation,  and  zero  skin  friction.  Both  the 
boundary  layer  displacement  thickness  and  the  external 
velocity  are  now  unknown.   The  Mechul  function  method  is 
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used  to  solve  the  flow  under  these  conditions  by  writing  the 
edge  boundary  condition  as  a  sum  of  the  inviscid  velocity 
distribution  and  the  perturbation  velocity  due  to  viscous 
effects.  The  perturbation  velocity,  Sue(x) ,  is  determined 
from  the  Hilbert  integral  [Refs.  12,27] 

<5ue(x)  =  -J   4-  (u  6*)  ~-  (6.4) 

ev  '         it  ;   da   e    x-o  v    ' 

a 
where : 

d(u  <S*) 
g  e is  the  blowing  velocity. 

The  solution  method  follows  the  direct  boundary  layer  method 
with  several  exceptions.  The  Falknar-Skan  transformation, 
with  its  constant  boundary  layer  thickness,  is  no  longer 
applicable;  nor  is  using  ue  as  a  reference  velocity  since  ue 
is  now  unknown.   Instead, 


n.  =  \fe   y  (6.5a) 


VX 


f(x,n)  = —  <Mx,y)  (6.5b) 

/u  VX 
o 


where  the  reference  velocity  is  now  taken  as  u0,  the  free 
stream  velocity.  The  solution  is  again  a  downstream 
marching  technique.  The  FLARE  (Flugge-Lotz  and  Reyhner) 
approximation  is  used  to  continue  the  integration  through 
regions  of  backflow.   In  these  regions,  where  the  velocity 
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in  relation  to  the  forward  velocities  is  assumed  small,  the 
streamwise  convection  term  u8u/8x  is  set  equal  to  zero. 
3.   Transition  Model 

For  the  IBL  technique  to  be  successful,  the 
displacement  thickness  must  be  accurately  determined.  This 
is  dependent  on  an  accurate  solution  of  the  laminar  and 
turbulent  flow  equations,  the  transition  region  between 
them,  and  when  applicable,  separated  flows. 

The  values  of  flow  parameters  associated  with 
laminar  and  turbulent  flow  differ  greatly:  the  boundary 
layer  thickness,  momentum  thickness,  skin  friction,  velocity 
profile,  and  drag  are  all  influenced  by  increased 
turbulence.  For  a  code  to  accurately  model  the  boundary 
layer  flow,  both  regimes  must  be  modeled  as  well  as  the 
transition  between  the  two  and  separated  flow.  The 
influence  of  the  transition  location  and  length  of 
transition  on  the  accuracy  of  the  solution,  especially  for 
low  Reynolds  flows,  has  been  demonstrated  by  Reference  12. 

The  Cebeci  IBL  code  uses  Michel's  criterion  to  predict 
onset  of  transition.  Michel  begins  transition  when  the 
local  Reynolds  number  based  on  the  momentum  thickness,  RQ/ 
is  related  to  the  length  Reynolds  number  by  the  empirical 
equation 


22.400     0.46 
Re  =  1.174(1  +  — p" —  )  R  (6.6) 
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where : 

Re  =  ueelv  ;  and 
Rx  =  uex/v  . 

The  transition  model  is  also  highly  empirical  and  is 
as  previously  given  in  Equation  (3.20). 

As  implied  by  the  above  discussion,  turbulence  onset 
and  its  generation  mechanisms  are  not  thoroughly  understood. 
Nonetheless,  it  is  evident  that  numerical  methods  need  to 
include  transition  capabilities  to  begin  to  accurately  model 
the  boundary  layer  flow. 

B.   INVISCID  OUTER  FLOW 

Hess  and  Smith  developed  a  panel  method  where  the  flow 
is  represented  by  a  series  of  vortices  and  sources.  [Ref. 
30]  They  assume  the  vortex  strength  to  be  constant  and 
distributed  over  the  surface  such  that  the  correct 
circulation  results.  The  velocity  field  is  then  modeled 
through  a  source  distribution  that  forces  the  velocity  to  be 
everywhere  tangent  to  the  surface  (the  no  penetration 
condition) .  This  method  is  simplified  by  defining  nodes  on 
the  airfoil  surface  and  connecting  them  with  straight  line 
panels.  Obviously,  greater  accuracy  is  achieved  by 
increasing  the  number  of  panels  on  the  airfoil.   [Ref.  12] 

The  Kutta  condition  must  also  be  satisfied.  The  Kutta 
condition  insures  a  unique  solution  by  imposing  zero  loading 
in  the  region  of  the  trailing  edge.    The  three   basic 
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principles  of  the  Kutta  condition  are  that  the  circulation 
about  the  airfoil  is  such  that  the  flow  leaves  the  trailing 
edge  smoothly,  the  trailing  edge  is  a  stagnation  point  if 
the  trailing  edge  is  finite,  and  the  upper  and  lower 
trailing  edge  velocities  are  finite  and  equal  if  the 
trailing  edge  is  cusped.  The  vortex  strength  is  determined 
from  the  Kutta  condition  and  the  source  strength  can  then  be 
calculated  from  the  vortex  strength.   [Ref.  30] 

C.   STRONG  INTERACTION  MODEL 

The  inner  and  outer  flow  influence  each  other  and  thus 
cannot  be  solved  separately  if  viscosity  effects  on  pressure 
are  large.  The  strong  interaction  model  couples  the 
boundary  layer  and  the  external  viscous  flow  by  allowing 
both  the  displacement  thickness  and  the  pressure  (which  is  a 
function  of  external  velocity)  to  be  unknown.  An  iterative 
simultaneous  solution  is  then  achieved  by  alternating 
between  the  viscous  and  inviscid  flow  equations  until 
convergence  is  achieved. 

The  solution  method  is  based  on  conformal  mapping  and 
Fourier  analysis  techniques  [Ref.  31].  The  airfoil  is 
mapped  onto  a  circle  through  a  series  of  conformal  mappings 
and  application  of  the  fast  Fourier-transform  algorithm.  It 
then  is  mapped  onto  another  circle  that  includes  the 
boundary  layer  and  so  models  the  inviscid  portion  of  the 
flow  as  that  over  an  airfoil  whose  boundaries  have  been 
displaced  by  the  viscous  boundary  layer.   At  the  surface  of 
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the  circle,  the  no  penetration  condition  does  not  apply.  To 
account  for  this,  a  nonzero  normal  velocity  (blowing 
velocity)  is  prescribed.   [Ref.  29] 

The  boundary  layer  is  solved  from  the  surface  to  the 
outer  boundary.  The  outer  boundary  conditions  are  defined 
by  the  interaction  law 


K 
(Ufe)i,k  »  (u  )i,k-l  +  l     Cik[[Ue6*)k,K-l(Ue6*)k/K-l    (6.7) 

k=l 


The  solution  of  this  is  an  approximation  that  requires 
several  sweeps  over  the  upper  and  lower  surfaces  to  achieve 
convergence. 

Convergence  between  the  boundary  of  the  two  methods 
is  checked  and  the  procedure  is  repeated  by  updating  the 
product  of  the  external  velocity  and  displacement  thickness 
until  desired  convergence  is  achieved. 

D.   USE  OF  CEBECI'S  IBL  CODE 

Reference  29  details  the  use  of  the  Cebeci  IBL  code. 
The  code  is  submitted  to  the  Naval  Postgraduate  School  IBM 
mainframe  via  Job  Control  Cards  where  account  information, 
running  time,  and  size  commands  are  set. 

Output  information  includes  input  data,  coefficients  of 
lift  and  drag,  airfoil  coordinates,  shear  stress,  skin 
friction,  displacement  and  momentum  thickness,  and  velocity 
profile  information.  The  inputs  to  the  program  are  located 
in  data  cards  and  are  defined  in  Table  6.1.   See  Reference 
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29  for  a  detailed  program  description.   An  example  of  input 
parameters  used  is  given  in  Table  6.2 . 
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TABLE  6.1 
INPUT  PARAMETERS  FOR  THE  CEBECI  IBL  CODE 


INPUT 
IPT(l) 

IPT(2) 

IRSTRT 

IGLMAX 

I PRINT 

RL 

XTRL 

XTRU 

ALPO 

GTR 

MPTS 

XP 

YP 


DESCRIPTION 

Flag  for  specification  of  lower  surface 
transition 

Flag  for  specification  of  upper  surface 
transition 

Flag  for  restarting  solution 

Number  of  sweeps 

Flag  for  output  printed 

Reynolds  number  in  millions 

Lower  surface  specified  transition 

Upper  surface  specified  transition 

Angle  of  attack 

Empirical    constant    for    length    of 
transition 

Number  of  coordinates 

X  airfoil  coordinates 

Y  airfoil  coordinates 
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TABLE  6.2 
VALUES  USED  FOR  THE  CEBECI  IBL  CODE 


INPUT 

IPT(l) 

IPT(2) 

IRSTRT 

IGLMAX 

I PRINT 

RL 

XTRL 

XTRU 

ALPO 

GTR 

MPTS 

XP 

YP 


VALUE 

4 

1 

0 

16 

2 

1-15 

0 

.005 

0-12 

1200 

101 

X  airfoil  coordinates 

Y  airfoil  coordinates 
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VII.   PRESENTATION  OF  COMPUTATIONS 

Simulations  of  steady  flow  at  Reynolds  numbers  ranging 
from  1  million  to  15  million  at  zero  degrees  angle  of 
attack,  angle  of  attack  studies  for  Reynolds  numbers  of  1.5 
million  and  6  million  all  over  a  NACA  0012  airfoil  are 
presented  in  this  chapter. 

Four  aerodynamic  factors  are  investigated:  coefficient 
of  lift,  coefficient  of  pressure,  velocity  profiles  at 
specified  locations  along  the  chord  of  the  airfoil,  and  skin 
friction  along  the  airfoil  chord.  In  the  steady  flow  cases 
results  from  the  Sankar  N-S  code  are  compared  with  results 
from  the  Cebeci  IBL  code.  Numerous  studies  have  documented 
the  validity  of  the  Cebeci  code  in  steady  flow  [Refs.  2,12, 
28,30],  however,  the  Cebeci  IBL  code  only  considers 
incompressible  flow.  Therefore,  the  Sankar  N-S  code  was 
limited  to  the  upper  regions  of  what  is  normally  considered 
incompressible  flow,  .3  Mach.  Unless  otherwise  indicated, 
all  cases  are  run  at  .3  Mach. 

The  coefficient  of  lift  as  a  function  of  angle  of  attack 
is  plotted  in  Figure  7.1  for  a  NACA  0012  airfoil  at  a 
Reynolds  number  of  6  million.  Figure  7.2  shows  lift 
coefficients  versus  angle  of  attack  for  a  Reynolds  number  of 
1.5  million. 
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Next,  the  correlation  between  the  suction  pressure 
coefficients  for  various  Reynolds  numbers  and  angles  of 
attack  were  investigated.  At  zero  degrees  angle  of  attack 
the  coefficient  of  pressure,  Cp  of  the  upper  surface  is  the 
same  for  the  five  Reynolds  numbers  presented:  1,  3,  6,  10, 
and  15  million,  using  both  Cebeci' s  IBL  code  and  Sankar's 
N-S  code.  Figure  7.3  plots  this.  Coefficients  of  pressure 
at  angles  of  attack  of  0,  4,  8,  and  12  degrees  at  a  Reynolds 
number  of  6  million  and  angles  of  attack  of  0,  2,  4,  and  6 
degrees  for  a  Reynolds  number  of  1.5  million  and  a  Mach  of 
.12  are  shown  in  Figures  7.4  and  7.5  respectively.  The 
conditions  of  Figure  7.5;  1.5  million  Reynolds  number,  .12 
Mach,  and  0,  2,  4,  and  6  degrees  angle  of  attack,  were 
chosen  to  verify  steady  state  conditions  presented  in  Tang 
[Ref.  23].  The  pressure  coefficients  were  compared, 
however,  no  other  results  were  presented  by  Tang. 

Skin  friction  and  velocity  profile  information  was  then 
sought.  The  verification  studies  done  by  References  3,  7, 
14,  and  15  did  not  address  either.  Therefore,  the  only 
comparisons  made  were  with  the  Cebeci  IBL  code.  Figure  7.6 
shows  the  coefficient  of  skin  friction,  Cf,  as  a  function  of 
airfoil  chord.  The  Sankar  N-S  code  is  calculated  as 
turbulent  over  the  entire  airfoil  [Ref.  13].  The  Cebeci  IBL 
code  allows  specification  of  transition   from  laminar  to 
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turbulent  flow  or  allows  the  option  to  have  the  code  compute 
transition  [Ref.  30],  Figure  7.6  plots  Cf  for  the  NACA 
0012  airfoil  at  a  Reynolds  number  of  6  million  at  zero 
degrees  angle  of  attack  for  eight  cases,  three  using  the 
Cebeci  IBL  code  and  five  using  the  results  of  the  Sankar  N-S 
code: 

Cebeci  IBL  Code 

1)  turbulent  flow;  transition  at  .005c 

2)  computed  transition  at  .27c 

3)  laminar  flow;  transition  at  .70c 

Sankar  N-S  code 

4)  first  grid  point  at  .0001c 

5)  first  grid  point  at  .00005c 

6)  first  grid  point  at  .00002c 

7)  first  grid  point  at  .00001c 

8)  first  grid  point  at  .000005c 

The  C-grid  generates  41  points  in  the  direction, 
clustering  them  near  the  wall  and  stretching  them  further 
from  the  wall.  Specifying  the  first  grid  point's  location 
from  the  wall  determines  where  the  grid  clustering  begins. 

This  study  was  repeated  for  a  Reynolds  number  of  3 
million  at  zero  degrees  angle  of  attack  in  Figure  7.8.  Six 
cases  are  presented,  two  using  the  Cebeci  IBL  code  and  four 
using  the  results  of  the  Sankar  N-S  code: 
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Cebeci  IBL  code 

1)  turbulent  flow;  transition  at  .005c 

2)  computed  transition  at  .44c 

Sankar  N-S  code 

3)  first  grid  point  at  .0001c 

4)  first  grid  point  at  .00005c 

5)  first  grid  point  at  .00001c 

6)  first  grid  point  at  .000005c 

Figure  7 . 8  presents  a  comparison  of  five  cases  for  a 
Reynolds  number  of  15  million  at  zero  degrees  angle  of 
attack: 

Cebeci  IBL  Code 

1)  turbulent  flow,  transition  at  .005c 

2)  specified  transition  at  .70c 

Sankar  N-S  code 

3)  first  grid  point  at  .0001c 

4)  first  grid  point  at  ,00005c 

5)  first  grid  point  at  .00002c 

Velocity  profiles  are  plotted  in  Figures  7.9  through 
7.17.  Figures  7.9  through  7.12  show  velocity  profiles  for 
the  conditions  of  Figure  7.6:  Reynolds  number  =  6  million, 
angle  of  attack  =  0,  N-S  Mach  =  .3.  Figure  7.9  is  the  IBL 
velocity  profile  computed  when  transition  is  specified  at 
.005c.  Figure  7.10  is  the  laminar  velocity  profile  resulting 
from  specifying  the  IBL  transition  at  .70c.  Figure  7.11  and 
7.12  are  the  velocity  profiles  computed  when  specifying  the 
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N-S  first  grid  point  from  the  wall  at  .0001c  and  .00002c 
respectively . 

Figures  7.13  through  7.15  are  the  velocity  profiles 
associated  with  the  conditions  of  Figure  7.8:  Reynolds 
number  =  15  million,  angle  of  attack  =  0,  and  N-S  Mach  = 
.30.  Figure  7.13  is  the  IBL  turbulent  velocity  profile, 
Figure  7 . 14  is  the  N-S  turbulent  velocity  profile  resulting 
from  specifying  the  first  grid  point  off  of  the  wall  as 
.0001c,  and  Figure  7.15  is  the  N-S  velocity  profile 
resulting  from  specifying  the  first  grid  point  off  of  the 
wall  as  .000005c. 

Figures  7.16  and  7.17  are  velocity  profiles  generated 
from  a  Reynolds  number  of  1  million  at  zero  degrees  angle  of 
attack  and  .30  Mach  number.  Figure  7.16  is  the  IBL  velocity 
profile  resulting  from  specifying  transition  at  .005c  and 
Figure  7.17  is  the  N-S  velocity  profile  resulting  from 
specifying  the  first  grid  point  at  .0001c  from  the  wall. 

The  results  of  Figures  7.6  through  7.8  are  replotted  in 
Figures  7.18  through  7.20.  However,  only  the  turbulent 
cases  and  small  grid  mesh  conditions  are  shown.  Figure  7.18 
plots  the  IBL  turbulent  skin  friction  at  zero  degrees  angle 
of  attack  for  Reynolds  numbers  of  3,  6,  and  15  million. 
Figure  7 .  19  plots  the  N-S  skin  friction  at  the  same  condi- 
tions as  Figure  7.19  for  the  results  of  the  case  with  the 
first  grid  point  specified  closest  to  the  wall.  Figure  7.20 
compares  the  IBL  and  N-S  results  from  Figures  7.18  and  7.19. 
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Next,  comparisons  between  the  IBL  and  N-S  skin  friction 
for  varying  angles  of  attack  were  investigated  for  the 
conditions  of  Figure  7.4:  Reynolds  number  =  6  million,  Mach 
=  .3,  and  Figure  7.5:  Reynolds  number  =  1.5  million,  Mach  = 
.12.  Figures  7.21  through  7.24  are  for  a  Reynolds  number  of 
6  million  and  a  Mach  of  .3.  Figure  7.21  plots  the  IBL  skin 
friction  profile  for  angles  of  attack  of  0,  4,  8,  and  12 
degrees.  Transition  is  specified  at  .005c.  Figure  7.22 
plots  the  N-S  skin  friction  for  the  same  angles  of  attack 
and  a  first  grid  point  off  the  of  the  wall  of  ,0001c. 
Figure  7.23  is  also  the  N-S  skin  friction,  but  the  first 
grid  point  off  of  the  wall  is  specified  at  .00001c.  Figure 
7.24  compares  the  IBL  and  N-S  skin  friction  profiles  for  12 
degrees  angle  of  attack. 

Figures  7.25  through  7.28  present  skin  friction  profiles 
for  a  Reynolds  number  of  1.5  million  and  a  Mach  of  .12. 
Figure  7.25  is  the  IBL  turbulent  skin  friction  (transition 
is  specified  at  ,005c).  Figures  7.26  and  7.27  are  the  N-S 
skin  friction  profiles  computed  when  the  first  grid  point  is 
specified  at  ,0001c  and  ,00005c,  respectively.  Figure  7.28 
compares  the  IBL  and  N-S  skin  friction  profiles  for  6 
degrees  angle  of  attack. 
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VIII.   RESULTS  AND  DISCUSSION 

The  Sankar  N-S  code  has  been  compared  to  the  Cebeci  IBL 
code.  Integrated  values  that  are  not  strongly  influenced  by- 
viscosity,  such  as  coefficients  of  lift  and  pressure, 
correlate  well,  as  shown  in  Figures  7.1  through  7.5,  and 
with  the  pressure  coefficients  presented  in  Tang  [Ref.  23]. 
However,  the  viscosity  influenced  boundary  layer  values, 
such  as  the  coefficients  of  skin  friction  and  velocity 
profiles,  are  very  sensitive  to  the  type  of  flow  and  the 
grid  size.  This  is  evident  in  Figures  7.6  through  7.8.  The 
IBL  laminar  skin  frictions  are  much  lower  than  the  turbulent 
values.  Also,  the  influence  of  the  grid  mesh  on  the  ability 
of  the  N-S  code  to  compute  turbulent  values  is  shown.  When 
the  first  grid  point  off  of  the  wall  is  chosen  as  .0001c, 
the  resulting  skin  frictions  appear  laminar  for  the  6 
million  and  15  million  cases  (Figures  7.6  and  7.8).  The 
lower  the  Reynolds  number,  the  less  sensitive  the 
computations  are  to  grid  size,  as  can  be  seen  by  comparing 
the  IBL  and  N-S  skin  friction  values  of  a  Reynolds  number  of 
3  million  (Figure  7.7).  In  all  cases,  specifying  the  first 
grid  point  closer  to  the  wall  positions  more  points  in  the 
boundary  layer.  For  Figure  7.6,  when  the  first  grid  point 
from  the  wall  is  specified  at  .0001c,  9  points  are  located 
in  the  boundary  layer  at  10%  chord,  16  points  at  50%  chord, 
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and  29  points  at  the  trailing  edge.  In  comparison,  when  the 
first  grid  point  from  the  wall  is  specified  at  ,000005c,  13 
points  are  located  in  the  boundary  layer  at  10%  chord,  19 
points  at  50%  chord,  and  29  points  at  the  trailing  edge  of 
the  airfoil. 

When  the  influence  on  skin  friction  of  subsequent  grid 
refinements  in  the  N-S  code  is  negligible,  the  N-S  skin 
frictions  are  substantially  greater  than  the  values  computed 
by  the  IBL  code  when  transition  is  specified  at  .005c.  See 
Figure  7.20.  It  is  of  interest  to  note,  however,  the  higher 
coefficients  of  skin  friction  evident  when  transition  is 
delayed  in  the  IBL  code  as  in  Figure  7.6. 

The  velocity  profiles,  Figures  7.9  through  7.17,  vary 
little  in  the  turbulent  cases,  regardless  of  grid  size  or 
Reynolds  number.  The  shape  of  the  profiles  varies  between 
the  IBL  and  N-S  cases,  but  both  exhibit  turbulent  velocity 
profiles.  There  is  a  discrepancy  between  the  skin  friction 
and  velocity  profile  results  computed  by  the  N-S  code  when 
the  first  point  off  of  the  wall  is  .0001c.  The  coefficients 
of  skin  friction  indicate  laminar  flow,  whereas  the  velocity 
profiles  are  turbulent.  A  laminar  velocity  profile,  Figure 
7.10,  is  calculated  using  the  IBL  code  (Reynolds  number  of  6 
million,  0  degrees  angle  of  attack) . 

The  influence  of  the  grid  mesh  is  again  evident  in  the 
angle  of  attack  studies,  Figures  7.21  through  7.28.  The 
coefficients  of  skin  friction  for  the  higher  Reynolds  number 
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of  6  million  seem  to  not  be  dependent  on  angle  of  attack  if 
the  grid  points  are  not  located  sufficiently  close  to  the 
wall.  This  is  not  evident  for  the  lower  Reynolds  number  of 
1  million.  The  N-S  skin  friction  values,  when  no  longer 
influenced  by  grid  size,  are  again  higher  than  the  fully 
turbulent  values  computed  by  the  Cebeci  IBL  code,  as  seen  in 
Figures  7.24  and  7.28. 

Varying  the  artificial  viscosity  in  the  Sankar  N-S  code 
did  not  vary  the  solution  of  the  skin  friction.  However, 
all  cases  were  run  with  the  first  grid  point  ,0001c  from  the 
wall,  so  no  conclusion  should  be  drawn.  If  the  number  of 
time  steps  was  exceedingly  large  (8000  steps)  for  zero 
degrees  angle  of  attack,  the  pressure  and  skin  friction 
profiles  began  to  indicate  separated  flow. 

Angle  of  attack  studies  using  the  Sankar  N-S  code 
required  a  large  number  of  time  steps  (approximately  7000) 
before  the  coefficient  of  lift  converged  to  realistic 
values,  regardless  of  the  first  grid  point  off  of  the  wall 
specified.  The  coefficient  of  friction  values  stabilized 
relatively  quickly,  within  2000  to  3000  time  steps.  Several 
computer  runs  at  each  angle  of  attack  and  Reynolds  number 
were  then  required  to  determine  when  specification  of  the 
first  grid  point  off  of  the  wall  no  longer  influenced  the 
skin  friction  values  computed.  At  this  point,  the  results 
were  assumed  to  be  converged. 
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The  IBL  code  took  less  than  16  iterations  to  converge,, 
regardless  of  the  Reynolds  number  chosen  or  the  angle  of 
attack.  Very  low  Reynolds  numbers  would  require  more 
iterations,  however. 
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IX.   CONCLUSIONS  AND  RECOMMENDATIONS 

Both  the  Sankar  Navier-Stokes  Code  and  the  Cebeci 
Interactive  Boundary  Layer  code  show  reasonable  results  for 
high  Reynolds  number,  incompressible  steady  flows  over  a 
NACA  0012  airfoil.  The  discrepancy  in  skin  frictions  should 
be  resolved  and  the  influence  of  dissipation  on  the  skin 
friction  investigated.  These  steady-state  results  should 
then  be  extended  to  other  airfoils  and  low  Reynolds  flows  to 
determine  the  effect  of  increased  viscosity  on  the  codes. 
The  effect  of  transition  on  the  velocity  profiles  and  the 
skin  friction  is  modeled  in  the  Cebeci  IBL  code;  however, 
the  lack  of  a  transition  model  and  a  smearing  function  in 
the  Sankar  N-S  code  limits  its  ability  to  model  most 
experimental  flows,  especially  at  low  Reynolds  numbers. 
Also  of  importance  is  the  wake  influence,  which  has  not  been 
addressed  in  this  report,  and  the  growth  of  the  boundary 
layer  along  the  airfoil.  Since  the  velocity  profile  results 
are  presented  non-dimensionally,  this  trend  is  not  evident. 
The  difference  in  profile  shapes  generated  from  the  N-S  code 
and  the  IBL  code  could  be  resolved  better  if  the  actual 
boundary  layer  profile,  rather  than  a  non-dimensional ized 
profile,  is  used. 

Neither  code  was  extensively  compared  to  experimental 
results  for  skin  frictions  or  velocity  profiles.   This  is, 


101 


of  course,  a  criterion  in  determining  the  accuracy  of  the 
codes.  The  cost  and  time  considerations  associated  with 
running  the  codes  indicate  that,  at  the  present,  for 
steady-state  runs,  the  IBL  code  is  more  effective.  The  N-S 
runs  were  submitted  to  the  NASA  X-MP  Cray,  where  waiting 
times  of  24  hours  before  the  program  was  executed  was 
common.  Execution  times  ran  between  10  and  40  minutes,  at  a 
cost  of  $1000  per  hour.  In  comparison,  the  IBL  program 
normally  was  completed  in  10  to  30  minutes  on  the  Naval 
Postgraduate  School  IBM  mainframe,  for  less  than  $50. 

The  N-S  code  is  a  very  effective  tool  for  calculating 
dynamic  stall  characteristics.  Its  capability  as  a  general 
flow  solver,  however,  is  limited  in  comparison  to  the  IBL 
technique  in  terms  of  cost  and  time  constraints. 
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